Approximating strongly correlated wave functions with correlator product states 
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We describe correlator product states, a class of numerically efficient many-body wave functions 
to describe strongly correlated wave functions in any dimension. Correlator product states in- 
troduce direct correlations between physical degrees of freedom in a simple way, yet provide the 
flexibility to describe a wide variety of systems. We show that many interesting wave functions 
can be mapped exactly onto correlator product states, including Laughlin's quantum Hall wave 
function, Kitaev's toric code states, and Huse and Elser's frustrated spin states. We also outline 
the relationship between correlator product states and other common families of variational wave 
functions such as matrix product states, tensor product states, and resonating valence bond states. 
Variational calculations for the Heisenberg and spinless Hubbard models demonstrate the promise 
of correlator product states for describing both two-dimensional and fermion correlations. Even in 
one-dimensional systems, correlator product states are competitive with matrix product states for 
a fixed number of variational parameters. 



How can one efficiently approximate an eigenstate 
of a strongly correlated quantum system? In one- 
dimensional systems, the density matrix renormaliza- 
tion group (DMRG) provides a powerful and system- 
atic numerical approach. 1,2 However, the accuracy of the 
DMRG in two or more dimensions is limited by the one- 
dimensional encoding of correlations in the matrix prod- 
uct states (MPS) that form the variational basis of the 
DMRGi 2 - Generalizations of MPS to higher dimensions — 
tensor network or tensor product states (TPS^^&i^ — 
have been introduced recently, but these engender con- 
siderable computational complexity (which does not arise 
with MPS). This has made it difficult to practically ex- 
tend the success and accuracy of the DMRG to higher 
dimensions. 

In this article we examine a different class of quan- 
tum states: correlator product states (CPS). Unlike MPS 
and TPS, which introduce auxiliary degrees of freedom 
to generate correlations between physical degrees of free- 
dom, CPS correlate the physical degrees of freedom ex- 
plicitly. The CPS form has been rediscovered many 
times,2iI2iii but the potential of CPS as an alternative 
to MPS/TPS for systematically approximating strongly 
correlated problems remains largely unexplored. Here 
we take up this possibility. CPS share many of the lo- 
cal properties of MPS/TPS but appear more suitable for 
practical calculations in more than one dimension as well 
as for fermion systems. To establish the potential of CPS, 
we analyze the relation between CPS and common fam- 
ilies of analytic and numerical trial wave functions. We 
then discuss the most important properties of CPS: they 
permit efficient evaluation of observables and efficient op- 
timization. Finally, we present variational Monte Carlo 
calculations for both spin and fermion systems. Our CPS 
results compare favorably with calculations using other 
variational wave functions that contain a similar number 
of variational parameters. 

Note: As this manuscript was completed we were in- 
formed of recent work by Isaev et al. on hierarchical 
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Figure 1: Nearest-neighbor 2-site and 2x2 plaquette CPS on a 
2D lattice. The CPS weight for a given quantum configuration 
\q\Qi ■ ■ ■ qi) is obtained by multiplying correlator coefficients 
together as in Eq. J2). 



mean-field theory— and by Mezzacapo et al. on entangled 
plaquette states^ as well as earlier work on string-bond 
states^ All these studies consider wave functions similar 
to CPS and share many of our own objectives. However, 
while our current efforts are related, especially to Ref. [H, 
we focus on aspects of CPS not addressed in these other 
works, such as the relationship with well-known analyti- 
cal and numerical wave functions, and we consider differ- 
ent physical problems, such as fermion simulations. Thus 
we regard our work as complementary rather than over- 
lapping. 



I. CORRELATOR PRODUCT STATES 

Consider a set of quantum degrees of freedom Q = 
{<Zi,92 ■ ■ on a lattice with L sites in one or more 
dimensions. Each qi might represent a spin S =1/2 de- 
gree of freedom, where q £ {|, J,}, or a fermion degree of 
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freedom, in which case q € {0, 1}. An arbitrary quantum 
state can be expanded over all configurations as 

\y) = J2* qiq2 - qL \qiq2---q L )- (l) 
{<?} 

A general quantum wave function requires an exponen- 
tial number of parameters — one for each configuration. 
One way to reduce the complexity of the problem is to im- 
pose some structure on the coefficients ^(Q). Correlator 
product states (CPS) are one example of this approach. 

CPS are obtained by associating variational degrees 
of freedom directly with correlations between groups of 
sites. For example, in the nearest-neighbor 2-site CPS, 
a correlator is associated with each neighboring pair of 
sites: 

m m 

where (ij) denotes nearest neighbors. The coefficients in 
Eq. ([T]) are given by products of correlator coefficients. 
For example, in a one-dimensional lattice, the amplitude 
of a configuration is 

*(Q) = C qiq2 C q2q3 C q3q4 . . . C qL - iqL . (3) 

Eq. ([2]) can be extended to higher dimensions simply by 
associating correlators with (overlapping) bonds on the 
lattice (Fig. []}. The nearest-neighbor 2-site CPS is an 
extremely simple CPS. Longer range correlations can be 
introduced by removing the nearest neighbor restriction 
on pair correlations or by including explicit correlations 
among more sites with correlators such as C qiq2q3 . It is 
clear that CPS provide a complete basis: in the limit of 
L-site correlators, the CPS amplitudes are precisely the 
coefficients of Eq. (fT]). 

When there is a global constraint on the total spin S 
or particle number N we can use projected CPS wave 
functions. For example, for fixed particle number, the 
A-projected nearest-neighbor 2-site CPS is 

M (ij) 

where Pjy ensures that g$ = A. Such projections 
do not introduce any complications in working with CPS 
and may be included in both deterministic and stochastic 
calculations without difficulty. 

It is sometimes useful to write the CPS in a different 
form. Each correlator element C qiqj can be viewed as 
the matrix element of a correlator operator C* y that is 
diagonal in the quantum basis {\qiqj)}: 

(q i q j \C»\qlq' j )=5 qigi> 8 qjql _C^. (5) 

The CPS wave function is obtained by acting a string of 
commuting correlator operators on a reference state |$). 
For example, a 2-site CPS may be written as 

|¥)=JJ£*|$>. (6) 

i>j 



When there are no constraints, the reference state is 
taken to be an equally weighted sum over all quantum 
configurations: |$) = £\i \q1q2 ■ ■ ■ otherwise, |$) 
is projected to satisfy the constraint. For example, if 
particle number is fixed, |<t> i\r) is an equally weighted sum 
over all quantum configurations with particle number A, 

\^>n) =J2^\qiq 2 ---qL)- (7) 

M 

Note that both projectors and correlators are diagonal 
operators in the Hilbert space and commute with one 
another: this means that the projection can be applied 
directly to the reference state and this simplifies numeri- 
cal algorithms using CPS. The operator representation is 
also useful when considering extensions to the CPS form 
such as alternative reference states. 



II. CONNECTION TO OTHER WAVE 
FUNCTIONS 

Many strongly correlated quantum states can be rep- 
resented exactly as correlator product states. CPS also 
have much in common with several classes of widely used 
variational wave functions: matrix product states, ten- 
sor product states, and resonating valence bond states. 
In this section, we discuss the connections between these 
wave functions. 



Huse-Elser wave functions 

In their study of frustrated spin systems, Huse and 
Elser constructed states in which the quantum ampli- 
tudes ^(Q) correspond to classical Boltzmann weights 
exp(— /3E[Q]/2) multiplied by a complex phased The 
weights are derived from an effective classical Hamilto- 
nian H cl . For example, in the case of pairwise correla- 
tions, H cl = EijHj with K] = K ijSi$l- The corre- 
sponding wave function can be represented as a 2-site 
CPS with C % i — exp(—phfj/2 + i(j)ij) where <pij assigns a 
complex phase to the pair ij. 

For the square and triangular Heisenberg lattices, Husc 
and Elser demonstrated that a very compact variational 
ground-state could be obtained with a semi-analytic 
three-parameter model for H cl (containing up to 3-site 
interactions) and an analytically determined phase. Al- 
though CPS can represent such highly constrained wave 
functions for symmetric systems, it can also serve as the 
foundation of a more general numerical method. By al- 
lowing correlators to vary freely and by considering hi- 
erarchies of larger correlated clusters, we can hope to 
construct rapidly converging approximations to arbitrary 
strongly correlated quantum states, as the DMRG does 
for one-dimensional quantum problems. 
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Laughlin wave function 

In 1983, Laughlin proposed a variational wave func- 
tion to explain the fractional quantum Hall effect.— The 
Laughlin wave function describes a strongly interacting 
system with topological order. Like the Huse and Elser 
wave functions, Laughlin's wave function can be associ- 
ated with the Boltzmann weights of an effective classical 
Hamiltonian and can be represented exactly as a corre- 
lator product state. 

The Laughlin quantum Hall state at filling fraction 
1 jm can be written in first quantization as 



of the toric code wave function can be generated by a 
CPS with plaquette correlators: 



ijkl _ I 1 ^ ^ l z SiS^S l z > 0, 

io if sisisisi < 0. 



(12) 



The exact representation of the toric code and Laugh- 
lin's wave function demonstrate the ability of CPS to 
describe systems with topological order. 



MPS and TPS 



*(ri,. 



■ riv) 



N 



(8) 



where z\ is the (complex) coordinate of particle A. (A 
Greek subscript indicates the coordinate of a particular 
electron. A Roman subscript indicates the coordinate of 
a lattice site.) Alternatively, the system can be mapped 
onto a discrete set of coordinates z± , . . . , Zl with an asso- 
ciated set of occupation numbers q\, . . . , Then Eq. © 
can be exactly expressed as a 2-site CPS in the occupa- 
tion number representation 



i*}=£n^n c 2 <93 ' p >i?i---^> 

{q} i i<3 
1 



Ci = 



c 2 = 



1 



1 



i (zi - z 3 y 



(9) 
(10) 
(11) 



The CPS wave function exactly reproduces the Laugh- 
lin wave function. It is, in some ways, more general than 
Eq. ([5]). The CPS form could be used to extend the 
Laughlin state beyond 2-site correlators while maintain- 
ing antisymmetry of the state, or to find a better varia- 
tional energy in open or disordered systems. 



Toric code 

Kitaev's toric code is another interesting quantum 
state with an exact CPS representation. Kitaev proposed 
the toric code as a model for topological quantum com- 
puting. The Hamiltonian is a sum of site and plaquette 
projectors on a square lattice with spins placed on the 
bonds. On a torus, the ground state of this Hamiltonian 
is 4- fold degenerate with a gap to all other excitations. 16 
It is an example of a quantum system with topological 
order. 

The ground state can be obtained from the zero- 
temperature Boltzmann weights of a classical Hamilto- 
nian H^ l mic — y^ n , hni- The sum is over all plaquettes 

Qj, and h 0i is a product of S z operators associated with 
the spins on the edges of the plaquette A The amplitudes 



Correlator product states are conceptually related to 
matrix and tensor product states. All of these wave func- 
tions can easily express entanglement between local de- 
grees of freedom. Nonetheless, CPS and MPS/TPS form 
different classes of quantum states and one is not a proper 
subset of the other. 

A matrix product state (MPS) is obtained by approx- 
imating the quantum amplitudes ^(Q) in Eq. |T]) as a 
product of matrices, one for each site on the lattice: 



HQ) 



w 



AOS 

12 12«3 



AlL 



(13) 



The "auxiliary" indices {«} are contracted in a one- 
dimensional pattern — a matrix product — and this gives 
rise to the low computational cost of working with MPS. 
However, the one-dimensional structure prevents MPS 
from efficiently describing correlations in higher dimen- 
sions. Tensor product states (TPS) extend MPS by ap- 
proximating the amplitudes ^{Q) by more general tensor 
contractions. Because of the more complicated contrac- 
tion pattern, TPS can in principle describe higher dimen- 
sional correlations^^ Unlike MPS, the TPS contraction 
cannot be evaluated efficiently in general. This leads to 
the high computational cost of working with TPS. 

To demonstrate the relationship between CPS and 
MPS/TPS, we consider a simple example of a nearest- 
neighbor 2-site CPS on a three-site lattice with periodic 
boundary conditions. This CPS amplitudes are 



\]/<?i92<?3 _ (jqiqi (jqiqz (jqzqi 



(14) 



Applying singular value decomposition to one of the cor- 
relators gives 



C qq 



£W 



(15) 



where we have absorbed the diagonal matrix Oi into Wf . 
With this decomposition, can be mapped to a MPS 
of auxiliary dimension 2: 



(W£U£)W£. 



(16) 
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This is equivalent to Eq. (13) with Afj = W?U]. The 
matrices of the resulting MPS (of dimension 2) have a 
restricted form. More complicated CPS (e.g., with 3- 
site correlators) map to MPS with larger auxiliary di- 
mension and more flexible forms for the matrices. (The 
dimension of the matrices grows exponentially with the 
range or number of sites in the correlator.) An arbitrary 
MPS cannot be mapped onto a CPS with less than the 
complete basis of L-site correlators. Conversely, a one- 
dimensional CPS with long-range correlators (such as the 
general 2-site CPS used in the Laughlin state) can only 
be represented by a MPS with an auxiliary dimension 
that spans the full Hilbert space. These arguments can 
be extended to higher dimensions and similar conclusions 
hold for the mappings between CPS and TPS. For a given 
number of variational degrees of freedom, only a subset 
of CPS can be exactly written as MPS/TPS and vice 
versa. 

While the correlators in the CPS have no auxiliary in- 
dices, they could be augmented by additional auxiliary 
indices. For example, string-bond states may be consid- 
ered one-site correlators with a pair of auxiliary indices.— 
n-site correlators can be generalized in an analogous way. 

The concept of an area law is sometimes used in the 
analysis of wave functions. If the amount of entangle- 
ment between a system and its environment scales with 
the area of the boundary between the two, the system 
is said to obey an area law. Arguments from quantum 
information theory suggest that wave functions that sat- 
isfy an area law can accurately describe systems (in any 
dimension) with a finite correlation length^ (Some crit- 
ical systems with long-range correlations also satisfy an 
area law, but others may violate the area law at zero 
temperature.) 

MPS wave functions satisfy a one-dimensional area law 
and have a finite correlation length. (Long-range correla- 
tions can be reproduced over a finite range, but they will 
eventually decay exponentially.) TPS satisfy area laws 
in two or more dimensions. CPS with local correlators 
like nearest neighbor pairs or plaquettes also satisfy an 
area law, making them promising candidates for systems 
with a finite correlation length. CPS with long-range cor- 
relators, such as those used in the Laughlin wave func- 
tion, are not constrained by an area law and can describe 
even more entanglement between system and environ- 
ment, obeying a volume law instead. 



RVB states 

Resonating valence bond (RVB) states are widely used 
in strongly correlated quantum problems i 18 i 19 A fermion 
RVB state can be written as a product of a Jastrow factor 
and a projected BCS wave function 

|*rvb) = e E « P N e E « A « a H |vac) (17) 

where Jy and Ay are commonly taken to be real. 



There is a close relationship between CPS and RVB 
states. At half-filling, the iV-projected 2-site CPS can 
be expressed in the form of Eq. (TTTj) . Consider a dimcr 
covering of the lattice. Let Ay = 1 for each pair ij that 
is connected by a dimer and Ay = otherwise. The 
corresponding projected BCS state is the CPS reference 
|$iv) defined earlier: 

P N e I^<X«J|vac) =P N J2\Q^2...q L ) = (18) 

{?} 

If the Jastrow factor is allowed to become complex, then 
the 2-site correlator C % i is fully parameterized as 

C ij = exp(J + J\hi + J{hj + J^hifij), (19) 

and the CPS and RVB wave functions are identical. 

Despite the existence of a mapping between the two 
wave functions, the emphasis of the CPS parameteri- 
zation is quite different from that of commonly studied 
RVB states. For fermion RVB wave functions where Jy is 
real, the Jastrow factor is positive and the nodes of the 
fermion wave function are those of the projected BCS 
state. In general, such a wave function cannot be exact. 
In contrast, the CPS wave function can modify the nodes 
of the reference wave function |$jv) through the complex 
Jastrow factor. By using higher order correlators, the 
CPS state can therefore become exact. While the most 
flexible RVB/ CPS form would combine a complex Jas- 
trow factor with an arbitrary projected BCS reference, 
there are computational advantages to the simpler CPS 
reference, including the possibility to efficiently evaluate 
observables without the use of a stochastic algorithm. 20 

III. COMPUTATIONAL COST OF CPS 

To be useful in practical calculations, a variational 
wave function must allow efficient evaluation of expecta- 
tion values and optimization of its parameters. 

This combination of properties in matrix product 
states is responsible for the success of the density matrix 
renormalization group. The expectation value of typi- 
cal observables can be evaluated exactly in a time which 
is polynomial in the size of the system. Likewise, the 
amplitude of a given configuration can also be evaluated 
exactly in polynomial time. As shown in Eq. (|13[) . the 
amplitude of a configuration is the trace of the product of 
L independent mxm matrices, where m is the dimension 
of the auxiliary indices {i} and L is the number of lattice 
sites. The cost for evaluating the amplitude is 0(m 3 L). 

Tensor product states generalize the structure of MPS 
to higher dimensions, but numerical efficiency is lost. In 
general, TPS amplitudes cannot be evaluated exactly in 
polynomial time! Additional renormalization procedures 
must be used while performing the contractions, which 
introduces an error that depends on the system size. For 
fermions, such errors can result in amplitudes or expec- 
tation values incompatible with a fermion wave function 



5 



as well as a variational energy below the fermion ground 
state, a so-called A-rcpresentability problem. As a re- 
sult, only certain classes of TPS are capable of efficient 
polynomial simulation. 

Like MPS, correlator product states allow efficient, ex- 
act evaluation of wave function amplitudes and expecta- 
tion values. For example, the amplitudes of a pair CPS 
are ^{Q) = J\i<j C qiqj . The amplitude is a simple prod- 
uct of numbers. This is true for any CPS, and thus the 
complexity is proportional only to the number of corre- 
lators in the ansatz. This is manifestly polynomial in the 
system size. In general, evaluation of the amplitude with 
n-site correlators will require O(L) multiplications if the 
correlators act locally — e.g, nearest neighbors, plaque- 
ttes, etc. — and 0(L n ) if there are no restrictions. 

This property allows efficient Monte Carlo sampling of 
expectation values. (Deterministic algorithms can also 
be used but will be presented elsewhere^) Moreover, 
constraints such as fixed particle number or total spin 
are easily handled within the Monte Carlo algorithm by 
limiting the Metropolis walk to states that satisfy these 
constraints. The expectation value of an operator is given 
by (A) = EqP(Q)MQ), where P{Q) = |vf(Q)| 2 and 



A(Q)=E(Q|i|Q')|^ 



(20) 



The sum over Q' extends over only those Q' for which 
(Q\A\Q') ^ 0. As long as A is sparse in the chosen 
basis, its expectation value can be evaluated efficiently. 
If |\&) is local (e.g., nearest-neighbor pair CPS), a further 
simplification occurs for operators such as a\a,j, alajdkOLi, 
or Si • Sj. For these operators, most of the factors in 
^(Q) and ^(Q') are identical and cancel from the ratio 
so that the time required to evaluate the expecation value 
is independent of the system size and depends only on the 
number of Monte Carlo samples. 

As with MPS and TPS, we can take advantage of the 
product structure of CPS when minimizing the varia- 
tional energy and use an efficient local optimization or 
"sweep" algorithm. The energy is minimized with re- 
spect to one of the correlators while the others are held 
fixed, then repeated for each of the correlators in turn un- 
til the energy has converged. This algorithm is described 
in more detail in the next section. 



IV. SPIN AND FERMION SIMULATIONS 

We have implemented a pilot variational Monte Carlo 
code to optimize general spin and fermion CPS wave 
functions. In Tables Q] and Ql] we present results for two 
models of interacting spins and fermions: (i) the 2D 
square Heisenberg model defined by the Hamiltonian 



H 



(ij) 



J J > 



(21) 



and (ii) a ID spinless Hubbard model at half filling. This 
model is defined by the Hamiltonian 



H = ^2-t(a\aj +a]cn) 

(ij) 



UriiUi 



(22) 



Each site can only be occupied or unoccupied, and the 
energy U is the cost of placing two fermions on neigh- 
boring sites. We studied periodic and open boundary 
conditions for both the Heisenberg and Hubbard models. 



Optimization method 

We optimize the correlators by minimizing the vari- 
ational energy with a sweep algorithm. At each step of 
each sweep, a target correlator is updated while the other 
correlators are fixed. Because the wave function is 
linear in the target correlator coefficients, the derivatives 
of 1^} with respect to these coefficients define a vector 
space for the optimization. For instance, if the target 
correlator has elements C M , then the vector space is gen- 
erated by the basis where 



l*M> 



30 



(23) 



Any vector in this space defines a CPS wave func- 
tion: x corresponds to the wave function |W(x)) = 
/sum^l^^}. 

It is convenient to work in a slightly different basis 
in which one vector Xo corresponds to the current value 
of the target correlator and the other vectors Xj are or- 
thogonal to xo (but not necessarily to each other). The 
updated target correlator will be a linear combination of 
the x a where a £ {0,i}. 

We construct the Hamiltonian H. a p and the overlap 
matrix S a p in this space: 



Ha^=(*(x ct )|i?|*(x /3 )) 

S af} =(*(x a )|*(x /3 )), 

where H is the model Hamiltonian of Eq. (|2ip or 
We then solve the generalized eigenvalue problem 

H ■ C = XS • C, 



(24) 
(25) 



(26) 



where C is a linear combination of the x Q . The eigen- 
vector with the lowest eigenvalue defines the optimal tar- 
get correlator coefficients & that give the lowest energy 
when all other correlators are fixed. We sweep over all 
of the correlators one at a time until the energy stops 
decreasing. 

This defines a general sweep algorithm for optimizing 
CPS. However to converge the sweeps when the Hamilto- 
nian and overlap matrix are constructed via Monte Carlo 
sampling it is very important to minimize the stochastic 
error. Nightingale and Melik-Alaverdian, 22 and Toulouse 
and Umrigar^ defined efficient estimators for variational 
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Table I: Variational Monte Carlo energies (in units of J) using CPS for the 2D S — 1/2 Heisenberg model, including percent 
errors (AE). CPS [2] denotes nearest-neighbor 2-site correlators and CPS [nxn] denotes plaquette correlators. The "exact" 6x6 
and 8x8 energies are obtained from a stochastic series expansion MC calculation using ALPS. 21 Unlike matrix product states, 
correlator product states maintain good accuracy as the width is increased. 



Lattice 


CPS [2] 


AE 


CPS [2x2] AE 


CPS [3x3] 


AE 


Exact 


Periodic Boundary Conditions 


4x4 


-11.057(1) 


1.5% 


-11.109(1) 1.1% 


-11.2202(2) 


0.1% 


-11.2285 


6x6 


-23.816(3) 


2.6% 


-24.052(2) 1.6% 


-24.313(2) 


0.5% 


-24.441(2) 


8x8 


-41.780(5) 


3.1% 


-42.338(4) 1.8% 


-42.711(3) 


0.9% 


-43.105(3) 


Open Boundary Conditions 


4x4 


-8.8960(5) 


3.2% 


-9.0574(4) 1.4% 


-9.1481(2) 


0.5% 


-9.1892 


6x6 


-20.811(1) 


4.2% 


-21.176(1) 2.5% 


-21.510(1) 


1.0% 


-21.727(2) 


8x8 


-37.846(3) 


4.5% 


-38.511(2) 2.8% 


-39.109(2) 


1.3% 


-39.616(2) 



Table II: Variational Monte Carlo energies (in units of t) for the L-site ID spinless Hubbard model with repulsion U using 
periodic and open boundary conditions, including percent errors (AE). CPS[n] denotes n-site correlators; DMRG[m] denotes 
a DMRG calculation with m renormalised states. Since CPS and DMRG calculations are not directly comparable in terms 
of complexity, the approximate number of degrees of freedom per site (d.o.f.) is listed in the bottom row. (The numbers are 
exact in the limit of an infinite lattice.) Encouragingly, CPS are competitive with MPS for a comparable number of variational 
parameters. Exact energies are from m=500 DMRG calculations. 



L 


U 


CPS [3] 


AE 


DMRG [3] 


AE 


CPS [4] 


AE 


DMRG [4] 


AE 


Exact 


Periodic Boundary Conditions 


12 





-7.052(1) 


5.5% 


-7.165 


4.0% 


-7.213(1) 


3.4% 


-7.313 


2.0% 


-7.464 


12 


4 


-2.692(2) 


4.1% 


-2.577 


8.2% 


-2.756(1) 


1.8% 


-2.725 


2.9% 


-2.807 


12 


8 


-1.461(1) 


1.1% 


-1.430 


3.1% 


-1.474(1) 


0.2% 


-1.462 


1.0% 


-1.477 


24 





-14.432(2) 


5.0% 


-14.608 


3.8% 


-14.714(2) 


3.2% 


-14.832 


2.4% 


-15.192 


24 


4 


-5.34(1) 


5.1% 


-5.340 


5.1% 


-5.482(1) 


2.6% 


-5.403 


4.0% 


-5.626 


24 


8 


-2.929(2) 


0.8% 


-2.860 


3.2% 


-2.931(1) 


0.7% 


-2.900 


1.8% 


-2.953 


36 





-21.82(1) 


4.6% 


-22.035 


3.6% 


-22.21(1) 


2.8% 


-22.421 


1.9% 


-22.860 


36 


4 


-7.93(3) 


6.0% 


-8.127 


3.7% 


-8.17(1) 


3.2% 


-8.173 


3.2% 


-8.440 


36 


8 


-4.390(2) 


0.9% 


-4.302 


2.9% 


-4.400(1) 


0.7% 


-4.355 


1.7% 


-4.430 










Open Bom 


■dary Conditions 










12 





-7.204(1) 


1.3% 


-7.185 


1.5% 


-7.274(1) 


0.3% 


-7.265 


0.4% 


-7.296 


12 


4 


-3.748(1) 


4.0% 


-3.787 


3.0% 


-3.887(1) 


0.5% 


-3.894 


0.3% 


-3.905 


12 


8 


-2.847(2) 


4.6% 


-2.920 


2.2% 


-2.971(1) 


0.4% 


-2.981 


0.1% 


-2.984 


24 





-14.593(1) 


2.2% 


-14.609 


2.1% 


-14.767(1) 


1.1% 


-14.838 


0.6% 


-14.926 


24 


4 


-6.32(1) 


7.8% 


-6.543 


4.5% 


-6.687(1) 


2.4% 


-6.803 


0.7% 


-6.851 


24 


8 


-4.287(2) 


6.6% 


-4.414 


3.8% 


-4.498(2) 


2.0% 


-4.576 


0.3% 


-4.590 


36 





-21.978(2) 


2.6% 


-22.035 


2.3% 


-22.260(2) 


1.3% 


-22.421 


0.6% 


-22.562 


36 


4 


-8.83(3) 


9.1% 


-9.323 


4.0% 


-9.36(1) 


3.6% 


-9.625 


0.9% 


-9.713 


36 


8 


-5.660(2) 


7.3% 


-5.873 


3.8% 


-5.934(3) 


2.8% 


-6.078 


0.4% 


-6.104 




d.o.f 


8 




18 




16 




32 







Monte Carlo optimization, and we have used these to con- 
struct Ti. and S. For numerical stability, it is important 
to monitor the change in the variational parameters and 
reject extremely large changes during a single iteration^ 
For CPS, this can be achieved by adding a dynamically 
adjusted diagonal shift to Ti that penalizes large changes 
away from C M : 5Hoo = 0, SHu > 0. Using this sweep 
algorithm, we find that the variational energy of the CPS 
converges (within statistical error) in less than 5 sweeps. 

To obtain the numbers in Tables Q] and [Til we ran the 
linear optimization routine for each system through 3 or 
4 sweeps, after which the energy stopped decreasing and 
instead fluctuated within a small range of values. We 
chose one wavefunction (set of correlators) from the final 



sweep and calculated the energy and variance reported in 
the tables using a larger number of Monte Carlo samples 
than we used during the optimization procedure. 

Results 

Table UJ shows the optimized energies obtained for the 
2D square Heisenberg model. This model tests the ability 
of CPS to describe two-dimensional correlations. When 
open boundary conditions are used, the system is not 
translation invariant and requires the kind of general pa- 
rameterization of the CPS emphasized here rather than 
the more restricted forms used by Huse and Elscr. 10 
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The nearest-neighbor 2-site CPS (CPS [2]) has only 
four variational parameters per site and gives errors in 
the range of 3-5% for open boundary conditions and 1- 
3% for periodic boundary conditions. The error is rapidly 
reduced by increasing the correlator size. For example, 
for the 8x8 periodic model, going from pair to 2x2 to 
3x3 plaquettes, the error goes from 3.1% to 1.8% to 
0.9%. The rapid convergence of the error with the corre- 
lator size is consistent with the results of Mezzacapo el 
al. for hardcore boson systems with periodic boundary 
conditions. 13 

As discussed earlier, CPS with local correlators like 
those used in Table [T] satisfy an area law, which allows 
them to accurately simulate systems with a finite corre- 
lation length. However, the 2D Hciscnbcrg model is gap- 
less with long-range correlations, so we expect the error 
to increase as the size of the lattice increases. Nonethe- 
less, the energetic error of the CPS wave function with a 
fixed correlator size grows quite slowly as the lattice size 
is increased. This is not true of MPS, in which the num- 
ber of variational parameters per site required to achieve 
a given accuracy grows rapidly with the width of a 2D 
system. 

We performed a series of DMRG calculations for the 
Heisenberg model on the lattices in Table [T] with a range 
of values of m using ALPS^i The variational objects in 
the DMRG are mxm matrices. For periodic boundary 
conditions, m 35, 250, and 750 are required for 1 per- 
cent accuracy on the 4x4, 6x6, and 8x8 lattices respec- 
tively. The latter calculation, which utilizes about 1.1 
million variational parameters per site (neglecting sym- 
metry and conservation of quantum numbers), is to be 
contrasted with the much more compact description us- 
ing the CPS with 3x3 correlators, which corresponds to 
just 512 parameters per site. 

The spinless ID Hubbard model with periodic bound- 
ary conditions has nontrivial fermion correlations and 
cannot be mapped onto a local spin model. Conse- 
quently, this model tests the ability of the CPS to capture 
fermion correlations. In Table [Hi we compare 3-site and 4- 
site nearest-neighbor CPS energies (CPS [3] and CPS [4]) 
with DMRG calculations for m = 3 and m = 4 renor- 
malised states. DMRG calculations were carried out us- 
ing ALPSj^i For open boundary conditions, the error 
in the CPS energy is smallest in the noninteracting sys- 
tem and largest for an intermediate interaction strength 
(U = 4). For periodic boundary conditions, the CPS [4] 
errors range from less than 1% for the U=8 case to ap- 
proximately 3% for the free fermion system — a difficult 
limit for a locally entangled state. The DMRG energies 
follow the same trends. 

To make a meaningful comparison with the DMRG 
results, we also show the approximate number of vari- 
ational degrees of freedom per site in each ansatz. A 
DMRG[to] wave function has 0(2m 2 L) degrees of free- 
dom (2 mxm matrices at each site) whereas the CPS[n] 
wave function has 0(2 n L) degrees of freedom (an n-site 
correlator at each site). As a result, the CPS [4] wave 



function has a similar complexity to the DMRG [3] state. 
Depending on the boundary conditions and the length of 
the lattice, the exact number of degrees of freedom may 
be less than this estimate. For instance, when L = 12 for 
an open chain, the DMRG [3] wave function has about 
14.7 parameters per site and the CPS[4] wave function 
has 12. Comparing the CPS and DMRG calculations 
with similar numbers of variational parameters, we see 
that the CPS energies are indeed very competitive, es- 
pecially for periodic boundary conditions, where a CPS 
includes direct correlations between the ends of the chain. 

Minimizing the CPS energy is a nonlinear optimiza- 
tion problem and the sweep algorithm may not converge 
to the global minimum of the variational energy. We 
have repeated the optimization for different initial wave 
functions to avoid local minima. The DMRG algorithm 
can also converge to a local minimum for m — 3 or 4. 
We repeated each of these DMRG calculation 100 times 
with the same input and reported the lowest energy ob- 
tained in Table [TTJ Although convergence to local min- 
ima is possible in both CPS and DMRG calculations, we 
believe the results reported in Tables Q] and |TT] indicate 
the competitive accuracy of CPS as a general variational 
method. 



V. CONCLUSION 

In this paper, we evaluated correlator product states 
as a route to describing strongly correlated wave func- 
tions in any dimension. Our preliminary numerical stud- 
ies indicate that CPS can capture both nontrivial fermion 
correlations and two-dimensional correlations. Together 
with the analysis showing the connections between CPS 
and many interesting quantum states, this supports the 
intriguing possibility that CPS are sufficiently flexible to 
systematically approximate general strongly correlated 
spin and fermion problems in two or more dimensions. 

Nonetheless, many questions remain to be answered. 
For example, how well do CPS reproduce correlation 
functions? While properties are harder to obtain ac- 
curately than energies in variational calculations, our 
view is that so long as successive CPS[n] calculations 
form a sufficiently rapidly convergent approximation to 
the quantum state of interest, then accurate approxima- 
tions to correlation functions can be constructed, as in 
the case of DMRG calculations. Detailed investigations 
of such questions and the analysis of more complex sys- 
tems like the full Hubbard model or the t-J model will 
require more sophisticated numerical treatments and al- 
ternative numerical techniques such as deterministic eval- 
uation methods. We are currently exploring these areas. 

We thank C.L. Henley for bringing the work of Huse 
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